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Abstract. We study the transition from the collisionlcss to the hydrodynamic regime 
in a two-component spin-polarized mixture of 40 K atoms by exciting its dipolar 
oscillation modes inside harmonic traps. The time evolution of the mixture is described 
by the Vlasov-Landau equations and numerically solved with a fully three-dimensional 
concurrent code. We observe a master/slave behaviour of the oscillation frequencies 
depending on the dipolar mode that is excited. Regardless of the initial conditions, 
the transition to hydrodynamics is found to shift to lower values of the collision rate 
as temperature decreases. 
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1. Introduction 

Several experiments on cooling two spin-polarized states of a Fermi gas have proven 
that this setup is an important tool for investigating the collisional properties of such 
a quantum system pQ. The JILA group has recently performed experiments in which 
the collisionality of a 40 K mixture is tuned by either varying the atomic density QjJ E] or 
the off-resonant value of the inter-species scattering length jSlllj- Also the temperature 
of the gas has been used to drive its collisionality and the importance of Pauli blocking 
at very low temperatures has thereby been demonstrated pQ. The role of collisions 
has also been investigated in dipolar collective modes of a boson-fermion mixture [H] at 
temperatures in both the quantum and classical regimes. The behaviour of the collision 
rate inferred from these experiments plays a crucial role in the development of strategies 
to improve the evaporative cooling techniques. In addition, the collisionality also affects 
other observables such as the damping and diffusion processes in the quantum mixture. 
In Ref. [I] an anisotropic expansion of a cigar-shaped two-component Fermi gas has 
been observed. In contrast to the bosonic case, the anisotropy of the expansion at very 
low temperatures cannot be taken as a clear signature of fermionic superfluidity [Hj and 
could be due to purely collisional effects. 

While ergodic approximations on semiclassical Boltzmann equations have been 
used to describe the evaporative cooling process [Zj, a microscopic understanding of 
the transition towards hydrodynamics £Q or of the anisotropic expansion in the case 
of attractive interactions [E] needs a complete numerical study of the dynamics of the 
quantum gas in position and momentum space [§]. In this paper we describe the strategy 
that we have used to solve numerically the Vlasov-Landau Equations (VLE) for an 
ultracold two-component fermion gas with a fully three-dimensional (3D) algorithm. 
The two fermionic fluids are treated by a particle-dynamics approach [TUj . accounting 
for mean-field interactions and for collisions between the two species. Since the Pauli 
principle causes a saturation of phase space at very low temperature, we have developed 
a locally adaptive importance- sampling technique which allows us to select the colliding 
particles in a way which is by several orders of magnitude faster than in standard 
Monte Carlo techniques. Furthermore, the control of the occupancy of the unitary cells 
of phase space and a suitable choice of computational parameters permit us to avoid 
Pauli inconsistencies during the dynamics. 

The two fermionic gases are confined inside different harmonic traps. Regardless of 
the initial perturbation exciting the dipolar modes, we observe that as temperature is 
lowered, even through most collisions become forbidden classically and by the Pauli 
principle, the few collisions that occur still suffice to drive the mixture from the 
collisionless to the hydrodynamic regime. This effect has already been predicted in 
our previous study [9 , where the angular degree of freedom was taken into account 
via an effective weight in a 2D numerical code. The limitations of the 2D approach are 
illustrated in the present work by comparing the oscillation frequencies and the damping 
rate of the dipolar oscillations evaluated within the 2D approach with those obtained 
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by means of the present 3D scheme. We also find that the choice of the initial state 
modifies the spectra of the dipolar modes. In particular, if the symmetry of the initial 
configuration is broken and one species is initially at rest, kicking by the other species 
will induce oscillations at the characteristic frequency of the latter. 

The paper is organized as follows. In Sec. El we introduce the physical model and 
describe in detail the numerical method developed to solve it. In Sec. El we study the 
dipolar oscillations in various collisional regimes and compare the results with those 
obtained in our previous work. Concluding remarks are given in Sec. EJ 

2. The model and its solution in 3D 

The two fermionic components in external potentials V^(r) are described by the 
distribution functions f^'(r,p,t) with j — 1 or 2. They obey the kinetic equations 

dtf {j) + - ■ V r f^ - V r U® • V p /^ = CitifW] (1) 
m 

where the mean-field effective potential is U^(r,t) = V^(r) + gn^(r,t) with j 
denoting the species different from j. Here we have set K — 1, g — 2ira/m r with 
a being the s-wave scattering length between two atoms of different species and m r 
the reduced mass, and n^\r,t) is the density obtained by integrating f^\v,p,t) over 
momentum. 

Collisions between atoms of the same spin can be neglected at low temperature, so 
that in Eq. Q the term Cu involves only collisions between particles that are polarized 
in two different Zeeman states. We have: 

C 12 [f^] = <r £ A„A tfUftofPf® - fMf?fTti j) ] (2) 

P2,P3,P4 

with /(?) ee fW(r,p,t), /W ee 1 - /Ci), /f ee /W(r,p is t), /f EE 1 - /f . V 
is the volume occupied by the gas and the factors A p and A e are the usual delta 
functions accounting for conservation of momentum and energy, with the energies given 
by pf/2m j + U^(r,t). 

The equilibrium state of the mixture is given by the stationary solution of Eq. 
i.e. by the local Fermi-Dirac distributions 

at given temperature T = l/ksP, where fij is the chemical potential ensuring the 
normalization condition J f^\r,p)d 3 r d 3 p/h 3 = Nj. The particle densities entering 
are to be determined self-consistently by integration over momenta. The details of 
this calculation have been previously given in Ref. [TT] . 

2.1. Numerical method 

The numerical solution of the coupled VLE in Eq. is carried out inside a finite box 
of size L x x L y x L z and discretized with a spacing Axj. We describe each species 
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in the box by means of a set of N c computational particles (cp) at phase-space points 
{ri,pi} which are spread according to the distribution /. In this way each fermion is 
represented by N q = N c /Nj > 1 "quarks". 

The initial selection of the computational particles is made by direct Monte Carlo 
sampling of the equilibrium distributions. Moreover, since collisions are very sensitive 
to the statistical nature of the colliding particles, we enforce the Pauli principle in each 
phase-space volume of size h 3 by resampling the computational particles that would 
otherwise exceed the allowed number of quarks N q in that cell. This further control is 
needed since we are dealing with a finite-size sample of Fermi-Dirac distributions. 

The solution of Eqs. proceeds in two steps: (i) the Lagrangian evolution of each 
distribution in phase space due to mean-field and external forces, and (ii) the occurrence 
of binary collisions among particles of the two species a la Boltzmann. Once the initial 
state is prepared, the first stage of the evolution is carried out by means of a second- 
order symplectic integrator O EH], and we update the position and velocity of each 
particle i according to the scheme 



Here F(r) is the mean- field force F(r) = -VU^(r) with j = 1 or 2. To compute 
the mean-field forces we first calculate the average density at a given cell by counting 
the number of particles inside the cell. The density is then numerically differentiated 
to obtain the forces at the grid points. The forces so calculated are finally linearly 
interpolated to the positions of the particles. 

It is worthwhile remembering that Liouvillean evolutions are incompressible, so 
that a distribution function satisfying the Pauli principle at a given time will satisfy the 
same constraint afterwards. This is an exact property whose numerical counterpart is 
checked during the course of the simulation. 

The second step of the solution involves binary collisions. At variance from the 
Lagrangian evolution these are tracked on a coarser mesh of spacing of the order of 
the de Broglie wavelength (a "Pauli cell"). In this way we introduce a temperature- 
dependent scale which imposes a length scale on which the particles are effectively 
indistinguishable. This can be seen as an extension of the standard Boltzmann approach. 
Furthermore, collisions in different Pauli cells can be handled independently due to the 
locality of Eq. ©. 

Our method allows us to count the number of collisions that occur step-by-step 
both classically and quantum-mechanically. In each Pauli cell we add up the classical 
probability of collisions between every pair of particles and store the pairs whose 
probability is greater than a locally adapted threshold. This reduces the running time 
and the storage needs of the code by several orders of magnitude, and at the same 
time the threshold is adjusted in order to guarantee a correct supply of pairs at each 
time step. The quantum features of the mixture are included in the last step of the 
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collision algorithm where the Pauli suppression in Eq. (J2J) is applied. For every classically 
colliding pair we calculate the blocking terms (1 — fV') by analyzing the occupancy of 
the momentum states in the Pauli cell. We have explicitly verified that this procedure 
yields the correct Pauli-suppressed fraction of collisions for a weakly interacting gas at 
least down to T ~ 0.2 Tp PH- 

We turn now to the discussion of the choice of computational parameters which 
allows us a reliable implementation of the Lagrangian evolution. During the dynamics 
two unphysical processes may occur: (i) the presence of an excess of quarks in saturated 
phase-space cells, and (ii) a broadening of the atomic cloud inside the trap. No cell of 
volume h 3 contains more than N q quarks for each species if the initial state has been 
accurately prepared. In this case, the inconsistent location of quarks takes place after 
we move the particles. Indeed, both problems (i) and (ii) are related to the calculation 
of the accelerations induced by mean-field forces. As already mentioned, these forces 
are constructed first on the simulation grid and then interpolated to the actual position 
of each particle. Therefore, an accurate description requires small spacings Axj and a 
large number of computational particles per cell. 

Aiming at understanding the impact on the results of the different choices for the 
values of N g and Lj/Axj we have studied the axial motion of two interacting Fermi 
clouds for several values of these parameters. Hereafter we shall consider magnetically 
trapped 40 K atoms that are equally shared in two different Zeeman states (mj? = 9/2 and 
mp = 7/2) and confined in isotropic harmonic traps with slightly different frequencies 
(ug/2 = 2tt x 19.8 s -1 and = 2tt x 17.46 s -1 ). For this test case we take the total 
number of particles N = 200, the mutual scattering length a = 3 x 10 4 Bohr radii 
(a strongly collisional regime), and T = 0.3 Tp. Both clouds are suddenly displaced 
from their equilibrium positions along the axial direction and start oscillating. We then 
analyze the evolution of the mixture for a total time tfi n ~ 50aC, 2 (~ 8 complete 
oscillations) and from the value of the centre-of-mass position Zj(t) we extract the 
oscillation frequencies Uj of the two components and the average damping rate 7. The 
details of this procedure will be given in Sec. El In addition we also compute the axial 
spreads of the density profiles Dj = (Azj) tfin / (Azj) to and the quantum collision rate 
T q by direct counting of the total collisions per particle, and monitor the number M 
of computational particles located in forbidden phase-space cells. The latter shall be 
called Pauli inconsistencies. For the system parameters here considered, we find that 
the oscillation frequencies do not vary considerably as N q and Lj/Axj are varied. The 
dependence of the other observables is summarized in Table ^ 

The inconsistencies disappear by increasing the number of grid points, while a large 
number of computational particles per cell is required to prevent statistical noise from 
producing spurious diffusion of the particles. The average damping rate 7 is sensitive 
to the spreading of the clouds, but seems to be less affected by the presence of few 
inconsistencies. On the other hand the collisional rate T q depends more strongly on 
both factors. 

In summary, this test case has given us enough information to select appropriate 
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Table 1. The maximum relative value M = (Af/N c ) max of particles in excess per time 
step, the maximum spread D = max{_Dg/ 2! -D7/2} of the clouds in the axial direction, 
the collision rate T q (in units of 1/s), and the average damping rate 7 (in units of 
1/s) as functions of the number of quarks N q and of grid points Lj/Axj. The data 
correspond to the two-component Fermi gas at T = 0.3Tf with a total number of 
particles N = 200 and an inter-species scattering length a = 3 x 10 4 ao (cto being the 
Bohr radius). 
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parameters for reliable calculations. In the following numerical experiments we shall 
choose the most suitable values of N q and Lj/Axj in order to ensure Dj ~ 1 and a 
negligible A/", while keeping the computing time within reasonable limits. 

3. Dipolar oscillations 

The quantum features of the mixture are most conveniently displayed by its collective 
modes, which in turn can be investigated by analyzing the dynamical response of 
the system under distortions of the equilibrium state. This can be done quite 
straightforwardly in actual experiments where the cloverleaf setup is used []] 121 El E] : a 
bias magnetic field displaces the centre of the trap, while a Stern-Gerlach separation of 
the two Zeeman states could be obtained by adding a magnetic field gradient [H [TH] . 

In the following we analyze the dipolar modes of the mixture at varying scattering 
length, as can be done e.g. by exploiting Feshbach resonances 0, and focus on two 
specific types of initial conditions: one in which both components start oscillating with 
the same amplitude (in-phase oscillations) and the other in which only the mp = 9/2 
species is initially displaced (kicked oscillations). 

3.1. In-phase oscillations 

The equilibrium fermionic density profiles for 2 x 10 4 particles are first created with 
the traps centred at r — (0, 0, 1) ay lo . At t = we let the confinements return to the 
position r = and hence the atoms start oscillating. In the collisionless regime the 
natural oscillation frequencies of the two species are different and correspond to trap 
frequencies renormalized by mean-field interactions, whereas the hydrodynamic regime 
is characterized by a common oscillation frequency of the two clouds, whose motions 
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Figure 1. Transition from the collisionless to the collisional regime in a 40 K mixture 
with -/V7/2 = -/V9/2 = 10 4 atoms, (a) Oscillation frequencies (in units of l/s) as a 
function of the collision rate T q (in units of l/s). The circles correspond to T = 0.3 Tp 
while the squares to T = 0.6 Tp; filled symbols mark the mj? = 9/2 species and the 
empty ones the rap — 7/2 species. The dash-dotted lines indicate the location of 
topi pt while the horizontal dashed lines indicate the bare trap frequencies, (b) Average 
damping rate 7 (in units of l/s) as a function of T q . 



have been locked by frequent collisions between their particles. In this limit c lassical 
kinetic theory predicts that the common frequency is ujhd = y(^|/2 + ^7/2 which 
in our case gives ujhd — 117.3 s _1 . This is in agreement with our numerical results as 
is shown in Fig. ^a), where the in-phase dipolar frequencies are plotted as functions of 
the collision rate T q . 

In evaluating the oscillation frequencies and damping rates we have assumed that 
the dynamics is well described by an exponentially damped harmonic motion and 
fitted each centre-of-mass position Zj(t) with a function Ajcos(ujjt + <f>j) exp(— jjt). 
We calculate the average damping rate as 7 = (79/2 + 77/2 )/2 and estimate the 
uncertainty of the result by analyzing the dependence of the fit parameters on the 
fitting range. This uncertainty is reported in the figures as error bars. The intermediate 
regime is characterized by a strong damping and larger uncertainties in the single-mode 
frequencies, which destroy the coherence of the oscillations. 

As an illustration of the various dynamical regimes we show in Fig. |2 the evolution 
of the centre of mass of the two species for three values of T q at T = 0.3 Tp. For 
27iT q ujJ 1 <C 1 the two clouds oscillate independently without appreciable damping, 
whereas as 2nT g ujJ 1 approaches unity the two centres of mass start moving together but 
still incoherently and therefore their oscillations are damped. At higher collisionality 
(27rr g cj7 1 > 1) the two species are locked and oscillate at the same frequency without 
damping. 

From Figures ^a) and^b) it is clear that the transition towards the hydrodynamic 
regime occurs for lower values of T q at T = 0.3 Tp than at T = 0.6 Tp. This phenomenon 



Transition to hydrodynamics in colliding fermion clouds 



8 




^9/2 t 



Figure 2. Ccntrc-of-mass position Zj(t) (in units of the oscillator length aho = 
as functions of u>g/2t for a mixture of 40 K atoms at T = Q.3Tp. From 
top to bottom, the indicated values of the collisional rate correspond to a = 150ao 
(collisionless regime), a = 5.0 x 10 3 ao (intermediate regime), and a — 1.5 x 10 4 a 
(hydrodynamic regime), in units of the Bohr radius ao- 

can be attributed to collisions which involve particles in a narrower energy range around 
the Fermi level as the temperature decreases jHj- However, the results presented here are 
only in qualitative agreement with those obtained in Refs. In fact, the 2D approach 
overestimates the collision rate and the damping rate. This is shown in Fig. El where 
the oscillation frequencies are compared with those obtained by using the 2D code. 

3.2. Kicked oscillations 

We next excite the dipolar mode with a different initial configuration, by taking the trap 
centre at (0, 0, 0) for the mp = 7/2 component and at (0, 0, 1) a^ for the mp = 9/2 
component. The confinement for the mp = 9/2 component is then shifted to z = 
and only its atoms start oscillating. If the two clouds did not interact, the m F = 9/2 
component would oscillate at its bare frequency and the centre-of-mass of the other 
species would not move. Collisionality drives both clouds into oscillation. 

In Fig. |IJa) we show the axial centre-of-mass position for the two clouds as a 
function of time for four values of the collision rate T q . At variance from what was 
found in the in-phase case, a single-mode model cannot describe the oscillations since 
even in the collisionless region the mp = 7/2 motion has a strong frequency component 
induced by the mp = 9/2 motion. This can be seen from the corresponding Fourier 
transform of the centre-of-mass positions in the various regimes, which are shown 

in Fig.Efb). In the case of small collisionality (the two top curves in Fig.0} the = 9/2 
component acts as a drive on the mp = 7/2 component but is not strongly affected by 
collisions. 
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Figure 3. Comparison of the oscillation frequencies ojj (in units of 1/s) obtained with 
the present 3D approach (squares) with those of an effective 2D system (circles) for 
a 40 K mixture with N g / 2 = N 7 / 2 — 100 particles at T — 0.3 If . The empty symbols 
correspond to the mp — 7/2 species and the filled ones to the mp = 9/2 one. 

The corresponding Zj(t) trajectories are drawn in Fig. ^a). For low collisionality 
(top and second row in Fig. a beating between the two peaks in the mp = 7/2 
spectrum is visible in the z 7 / 2 (t) signal. At higher collisionality the two clouds get 
locked and oscillate at a single frequency, which is intermediate between the two trap 
frequencies (see also the spectrum in Figllfb))). 

4. Summary and concluding remarks 

We have studied the transition from the collisionless to the hydrodynamic regime in a 40 K 
mixture at temperatures where the effects of Pauli blocking are noticeable. The Vlasov- 
Landau equations of two interacting fermionic fluids are solved numerically by means 
of a fully three-dimensional particle-dynamics approach incorporating mean-field and 
collisional interactions. The implementation of a locally-adaptive importance sampling 
technique during the collisional step has allowed us to obtain a method which is very 
efficient at low temperature, for realistic number of particles and in a wide range of 
collisional strength. By exploiting this code we investigate the dipolar collective modes 
and analyze the oscillation frequencies and damping rates of the system for various 
initial configurations. We find that the transition to hydrodynamics is shifted to lower 
values of the collision rate as temperature decreases. 

As a future application, this procedure shall permit us to perform numerical 
experiments on cooling and expansion dynamics. These studies will help characterize 
the mixture in the intermediate regime where neither purely collisionless nor collisional 
approaches are appropriate. In addition the extension to highly interacting systems 
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Figure 4. Kicked oscillations of the mp — 9/2 component (full lines) and the 
m F = 7/2 component (dashed lines) of a 40 K mixture at T — 0.32V with N = 2 x 10 4 
atoms, (a) Centre-of-mass positions Zj(t) (in units of the oscillator length ato) as a 
function of w g / 2 from top to bottom the data correspond to a = 150 a (collisionless 
regime), a = 1.5 x 10 3 ao (intermediate regime), a — 5.0 x 10 3 clq (intermediate regime 
approaching collisional) , and a = 1.5 X 10 ao (hydrodynamic regime). The amplitude 
of z 7 / 2 (t) has been rescaled by a factor 10 and 5 in the top and second row, respectively, 
(b) Corresponding Fourier transform |£j(a>)| (in arbitrary units and log scale) as a 
function of the frequency u> (in units of w g / 2 )- The vertical dashed lines indicate the 
location of the bare trap frequencies. 



in the so-called unitary limit of collisions jTHI shall provide more physical insight into 
recent experiments on attractive mixtures of fermionic atoms |17j . Efforts to pursue 
these further studies are currently underway and will be reported elsewhere. 
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